;********************************************WI-CM-1-1-MR-vcbc.nc
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;load library for  wetbulb_stull
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/heat_stress.ncl"
;The *area_conserve_remap *and *area_conserve_remap_Wrap* have been *deprecated*.
;The* ESMF regrid* <http://www.ncl.ucar.edu/Applications/ESMF.shtml> library is the recommended replacement. It uses a superior methodology.
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
;********************************************
 begin

;set baseline
 dirPath = "~/source_data/"
 dirPath2 = dirPath+"/outdoor/"

 mo  = (/"ACCESS-CM2","AWI-CM-1-1-MR","BCC-CSM2-MR","CMCC-CM2-SR5","EC-Earth3","GFDL-ESM4", \
               "IITM-ESM","KIOST-ESM","MIROC6","MIROC-ES2L","MPI-ESM1-2-HR","MRI-ESM2-0"/); 12 models
 case = (/"ssp585" /)
 ncase = dimsizes(case);
 nmod = dimsizes(mo)

 syr = 1990; including historical data to identify 1C warming
 fyr = 2099;
 csyr = sprinti("%0.4i",syr);
 cfyr = sprinti("%0.4i",fyr);
 nyr = fyr - syr + 1;
 ndays = nyr*365

;read dimension for 1 deg data
  dimPath=dirPath+"/domain/"
  landfrac = "sftlf_fx_360x180_1deg_gn.nc"
  area = "areacella_fx_360x180_1deg_gn.nc"

  fdim  = addfile(dimPath+landfrac, "r")
  fdim2 = addfile(dimPath+area, "r")
  lat := fdim->lat
  lon := fdim->lon
  nlat = dimsizes(lat);
  nlon = dimsizes(lon);
  garea = fdim2->areacella  ;m^2
  lfrac = fdim->sftlf/100; %
  wgt := garea*1e-6*lfrac ; sqrkm
  wgt@_FillValue = default_fillvalue("float"); make sure FillValue is set correctly for wgt
  copy_VarCoords(garea, wgt)
  printVarSummary(wgt)
  print("global total land area: "+sum(wgt)+" sq km")



 ;===Use precalculated global warming amount relative to 1980-2009 (based on ERA5, HadCRUT5, BerkeleyEarth, NOAAGlobalTemp) to classify warming levels for each CMIP6 model
  ;see Data_global-warming-amount-ERA5-30yr-CMIP6_12models.ncl
  filename1 = dirPath+"Global_warming_amount_10yrma_"+syr+fyr+"_12models_obs30yr.csv"; 10-year running median
  filename2 = dirPath+"Global_warming_amount_30yrma_"+syr+fyr+"_12models_obs30yr.csv"; 20-year running median

;---Read in file as array of strings so we can parse each line
  lines  = asciiread(filename1,-1,"string")
  lines2 = asciiread(filename2,-1,"string")
  k = 1 ;number of header lines
  nlines = dimsizes(lines)-k   ; skip header lines
  delim = ","
;---Read fields (columns)
  models  =           str_get_field(lines(k:),1,delim)
  exp1   =           str_get_field(lines(k:),2,delim)
  central_10yr =     str_get_field(lines(k:),3,delim)
  period10yr =       str_get_field(lines(k:),5,delim)
 ;read median temperature anomaly from file 1
  Ta_anomaly10yr = tofloat(str_get_field(lines(k:),7,delim))
  ;return to 2d (nmod,nyr); nmod is not 15 but 12 now, be careful
  Ta_anomaly10yr@_FillValue=1e+20

  central_30yr =    str_get_field(lines2(k:),3,delim)
  period30yr =       str_get_field(lines2(k:),5,delim)
  Ta_anomaly30yr = tofloat(str_get_field(lines2(k:),7,delim))
  Ta_anomaly30yr@_FillValue=1e+20



 ;------choose whether 10-yr or 30-yr running period
  Ta_anomaly = Ta_anomaly10yr; either 10-yr or 30-yr running median/max
  periods = period10yr
  central_year = central_10yr

  
 ;--identify 0.8-4.5 warming amounts for detecting of WoE in the range of 1-4
 levels := fspan(0.8,4.5,38); 0.8 to 4.5 degC by 0.1 degC
 nlev = dimsizes(levels)


 DATA = True; 
 if (DATA) then

 Gday_sun_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday_diff_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday_diff1_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday_diff2_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday_zero_ens = new((/nmod,nlev,nlat,nlon/),"float")
 TWday_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gtwday_ens = new((/nmod,nlev,nlat,nlon/),"float")
 hc_ens = new((/nmod,nlev,nlat,nlon/),"float")
 P_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Taday_ens = new((/nmod,nlev,nlat,nlon/),"float")

 Gday0d_sun_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday0d_diff_ens = new((/nmod,nlev,nlat,nlon/),"float")
 Gday0d_zero_ens = new((/nmod,nlev,nlat,nlon/),"float")
 TWday35d_ens = new((/nmod,nlev,nlat,nlon/),"float")


 do mm=0,nmod-1

   print(" ")
   print("process CMIP6 model: "+mo(mm))


   infile = dirPath2+"/G-daynight-"+syr+fyr+"-"+mo(mm)+"-vcbc-outdoor.nc"; outdoor work

   fi = addfile(infile, "r")
   Gday_sun  := fi->Gday_mx_sun;  annual maxima of daytime average
   Gday_diff := fi->Gday_mx_diff
   Gday_zero := fi->Gday_mx_zero
   Gday_diff1 := fi->Gday_mx_diff_L
   Gday_diff2 := fi->Gday_mx_diff_U

   Gday0d_sun  := fi->Gday0d_sun
   Gday0d_diff := fi->Gday0d_diff
   Gday0d_zero := fi->Gday0d_zero
   TWday := fi->TWday_mx; (yrs:yrd,:,:)
   TWday35d := fi->TWday35d; (yrs:yrd,:,:)
   printVarSummary(TWday)
   printMinMax(TWday, 1)

   f3 = addfile(dirPath2+"/TWday-TQRUPZ-3h-"+syr+fyr+"-"+mo(mm)+"-vcbc.nc", "r")
   taC  := f3->Ta_TWday; (yrs:yrd,:,:,:);
   rh  := f3->RH_TWday; (yrs:yrd,:,:,:);
   wind2m := f3->U2_TWday; (yrs:yrd,:,:,:);
   pa  := f3->P_TWday; (yrs:yrd,:,:,:);
   taC_day = dim_avg_n_Wrap(taC, 1)
   rh_day = dim_avg_n_Wrap(rh, 1)
   P_day = dim_avg_n_Wrap(pa, 1)

   tdC := dewtemp_trh(taC + 273.15, rh) - 273.15        ; C
   TW3h  := wetbulb(pa/100, taC, tdC)
   copy_VarCoords(taC, TW3h)
   TWday = dim_avg_n_Wrap(TW3h, 1) ;night time values already set to NA
   printVarSummary(TWday)
   printMinMax(TWday, 1); same as above

   ;==convert TW to Gtw with Eq.10 or Eq.14==
   ;--use dynamic wind concurrent with TWday to convert it to Gtw 
   hc:= 14.1*wind2m^0.5 ;forced convection function from Fig.3 in Fiala 1999
   hc:= where(wind2m .le. 0.1, 3.3, hc); set natural convection hc=3.3 W m-2 K-1 for U<=0.1 m/s as in Fiala
   copy_VarCoords(taC, hc)
   hc_day = dim_avg_n_Wrap(hc, 1)

   ;baseline skin temperature for conversion of TW to Gtw
   Ts =35 ;default for TW and Gtw used by Sherwood and Huber (2010)
   esat_skin = satvpr_water_bolton(Ts, (/0,1/))  ;esat at skin temperature 35 degC

   cp=1010 ;specific heat capacity of air (Joules/kg/°C)
   lambda = latent_heat_water(Ts, (/0, 0/), 1, False)
   ;1) convert TW to Gtw by Eq.10
   ;esat_tw := satvpr_water_bolton(TW3h, (/0,1/))
   ;gamma0 = (pa*cp)/(0.622*lambda)
   ;Gtw3h := -hc*(Ts - TW3h) - hc/gamma0*(esat_skin - esat_tw)

   ;2) or by Eq. 14: Gtw = -hc(1+lambda/cp*delta)(Ts-Tw) = k(Ts-Tw)
   ;delta_q := (0.622/pa)*1000*(2508*exp(17.3*Ts /(Ts + 237.3))/(Ts + 237.3)^2) ; be consistent with nonlinear scale in Fig.3
   TW := TWday
   pa := P_day
   hc := hc_day
   delta_q := (0.622/pa)*1000*(2508*exp(17.3*(TW+Ts)/2 /((TW+Ts)/2 + 237.3))/((TW+Ts)/2 + 237.3)^2); more accurate when TW<<Ts
   k := -hc*(1+lambda*delta_q/cp);
   Gtwday := k*(Ts - TW) 
   copy_VarCoords(TWday, Gtwday)


   ;*concurrent Ta used to mask out cool conditions.
   f4 = addfile(dirPath2+"/Gday_zero-TQRUPZ-3h-"+syr+fyr+"-"+mo(mm)+"-vcbc-outdoor.nc", "r")
   ta3h  := f4->Ta_Gday_zero; 
   ta_zero_day = dim_avg_n_Wrap(ta3h, 1) ;daily mean Ta (night time values included)
 

   ;==option 2) 10-yr or 30-yr average 
   ;note Raymond et al. used 99.9th percentile and 30-year maximum, Pal & Elthair and Im et al. used 30-year maximum of daily maximum
   ;we can show 99.9th percentile when comparing with historical period 1980-2019
   ;but use 10-yr/30-yr maximum for each warming level because some models do not have 30-year median Ta for 4degC warming before 2099
 
   ; Set the window size for calc the average of a period corresponding to each warming level
   window_size = 10; 20; 30  


   ;==find the 10-yr mean corresponding to each warming level
   do ii=0,nlev-1 ;
     print("")
     print("process warming level: "+levels(ii))
 
     Ta_anomaly_m = Ta_anomaly(ind(models .eq. mo(mm))); select the specific model
     period_m = periods(ind(models .eq. mo(mm))); 10-yr/30-yr periods for each model
     central_year_m = central_year(ind(models .eq. mo(mm)))
 
     tDif = abs(levels(ii) - Ta_anomaly_m)
     ind1 = minind(tDif); closest warming amount
     yr1 = tointeger(central_year_m(ind1))
     ;matching within +-0.05 degC tolerance
     if (.not.(ismissing(ind1)) .and. tDif(ind1) .lt. 0.05) then
       print("closest matching period of "+mo(mm)+" for this level: "+period_m(ind1)+" with "+Ta_anomaly_m(ind1))
       print("at central year: "+yr1); use year index to pick the value when syr of Ta_anomaly is different from model data


      ystart = yr1 - window_size/2 +1
      yend = yr1 + window_size/2
      if (ystart.lt.syr .or. yend.gt.fyr) then
        continue
      end if
      print("start:end year for moving statistics:"+ystart+"-"+yend)

      Gday_sun_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday_sun({ystart:yend},:,:), 0))
      Gday_diff_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday_diff({ystart:yend},:,:), 0))
      Gday_diff1_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday_diff1({ystart:yend},:,:), 0))
      Gday_diff2_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday_diff2({ystart:yend},:,:), 0))
      Gday_zero_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday_zero({ystart:yend},:,:), 0))
      TWday_ens(mm,ii,:,:) = tofloat(dim_avg_n(TWday({ystart:yend},:,:), 0))
      Gtwday_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gtwday({ystart:yend},:,:), 0))
      hc_ens(mm,ii,:,:) = tofloat(dim_avg_n(hc_day({ystart:yend},:,:), 0))
      P_ens(mm,ii,:,:) = tofloat(dim_avg_n(P_day({ystart:yend},:,:), 0))
      Taday_ens(mm,ii,:,:) = tofloat(dim_avg_n(ta_zero_day({ystart:yend},:,:), 0))

      Gday0d_sun_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday0d_sun({ystart:yend},:,:), 0))
      Gday0d_diff_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday0d_diff({ystart:yend},:,:), 0))
      Gday0d_zero_ens(mm,ii,:,:) = tofloat(dim_avg_n(Gday0d_zero({ystart:yend},:,:), 0))
      TWday35d_ens(mm,ii,:,:) = tofloat(dim_avg_n(TWday35d({ystart:yend},:,:), 0))

     end if
  end do

 end do; nmod



   ;======calculate ensemble median and confidence interval at each grid cell at each warming level
   Gday_sun_med = dim_median_n(Gday_sun_ens, 0)
   Gday_diff_med = dim_median_n(Gday_diff_ens, 0)
   Gday_diff1_med = dim_median_n(Gday_diff1_ens, 0)
   Gday_diff2_med = dim_median_n(Gday_diff2_ens, 0)
   Gday_zero_med = dim_median_n(Gday_zero_ens, 0)
   TWday_med = dim_median_n(TWday_ens, 0)
   Gtwday_med = dim_median_n(Gtwday_ens, 0)
   hc_med = dim_median_n(hc_ens, 0)
   P_med = dim_median_n(P_ens, 0)
   Taday_med = dim_median_n(Taday_ens, 0)
   print("mean ensemble median hc across 1-4 C warming used to convert TW to Gtw: "+avg(hc_med))
   ;print("mean ensemble median P across 1-4 C warming used to convert TW to Gtw: "+avg(P_med))


   Gday0d_sun_med = dim_median_n(Gday0d_sun_ens, 0)
   Gday0d_diff_med = dim_median_n(Gday0d_diff_ens, 0)
   Gday0d_zero_med = dim_median_n(Gday0d_zero_ens, 0)
   TWday35d_med = dim_median_n(TWday35d_ens, 0)

   ;=========calculate confidence interval at each grid cell 
    print("")
    ncrit = round(nmod/2.0, 3); round to integer 5 or 6 or 7 ;(half of models)
    print("nmod="+nmod+" | ncrit="+ncrit)

   ;--method 2) use the 25-75th percentile interval at each grid cell to account for inter-model uncertainty at specific locations e.g., megacities

   ;show 25–75th percentiles -> 50% confidence interval (often shown in other studies, Zhang et al. 2021, Matthews et al. 2017)
   ;?or show 16-84th percentiles -> 68% confidence interval (similar to 66% the Intergovernmental Panel on Climate Change’s ‘likely range’)

   nval := dim_num_n(.not.ismissing(Gday_sun_ens),0); Counts the number of True models at each warming level
   nval := where(nval.ge.ncrit, nval, default_fillvalue(typeof(nval))); only consider grid cells with >=6 models  
   nval@_FillValue = default_fillvalue(typeof(nval)) 
   printMinMax(nval,1)
   pctl_25 = round(0.25*nval,3)-1; 3rd low -> 3 models (0,1,2) at or below this value
   pctl_75 = round(0.75*nval,3)-1; 9th high -> 9 models (0-8) at or below this value
   ;printVarSummary(pctl_25)
   ;subtract 1 to get the correct model indices

   Gday_sun_sort = Gday_sun_ens
   Gday_diff_sort = Gday_diff_ens
   Gday_diff1_sort = Gday_diff1_ens
   Gday_diff2_sort = Gday_diff2_ens
   Gday_zero_sort = Gday_zero_ens
   TWday_sort = TWday_ens
   Gtwday_sort = Gtwday_ens
   Gday_sun_index = dim_pqsort_n(Gday_sun_sort, 2, 0); sorting models in increasing order at each cell
   Gday_diff_index = dim_pqsort_n(Gday_diff_sort, 2, 0)
   Gday_zero_index = dim_pqsort_n(Gday_zero_sort, 2, 0)
   TWday_index = dim_pqsort_n(TWday_sort, 2, 0)
   Gtwday_index = dim_pqsort_n(Gtwday_sort, 2, 0)

   Gday0d_sun_sort = Gday0d_sun_ens
   Gday0d_diff_sort = Gday0d_diff_ens
   Gday0d_zero_sort = Gday0d_zero_ens
   TWday35d_sort = TWday35d_ens
   Gday0d_sun_index = dim_pqsort_n(Gday0d_sun_sort, 2, 0); sorting models in increasing order at each cell
   Gday0d_diff_index = dim_pqsort_n(Gday0d_diff_sort, 2, 0)
   Gday0d_zero_index = dim_pqsort_n(Gday0d_zero_sort, 2, 0)
   TWday35d_index = dim_pqsort_n(TWday35d_sort, 2, 0)

   Gday_sun_low = new((/nlev,nlat,nlon/), typeof(Gday_sun))
   Gday_diff_low = new((/nlev,nlat,nlon/), typeof(Gday_sun))
   Gday_zero_low = new((/nlev,nlat,nlon/), typeof(Gday_sun))
   TWday_low = new((/nlev,nlat,nlon/), typeof(TWday))
   Gtwday_low = new((/nlev,nlat,nlon/), typeof(Gtwday))
   Gday_diff_high = new((/nlev,nlat,nlon/), typeof(Gday_sun))
   Gday_sun_high = new((/nlev,nlat,nlon/), typeof(Gday_sun))
   Gday_zero_high = new((/nlev,nlat,nlon/), typeof(Gday_sun))
   TWday_high = new((/nlev,nlat,nlon/), typeof(TWday))
   Gtwday_high = new((/nlev,nlat,nlon/), typeof(Gtwday))

   Gday0d_sun_low = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   Gday0d_diff_low = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   Gday0d_zero_low = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   TWday35d_low = new((/nlev,nlat,nlon/), typeof(TWday35d))
   Gday0d_diff_high = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   Gday0d_sun_high = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   Gday0d_zero_high = new((/nlev,nlat,nlon/), typeof(Gday0d_sun))
   TWday35d_high = new((/nlev,nlat,nlon/), typeof(TWday35d))

   do i=0,nlon-1
    do j=0,nlat-1 
     do ii=0,nlev-1
     ii_25 = pctl_25(ii,j,i)
     ii_75 = pctl_75(ii,j,i)
     if(ismissing(ii_25) .or. ismissing(ii_75)) then; 
        continue
     end if
     ;print(ii_25)

     Gday_sun_low(ii,j,i) = Gday_sun_sort(ii_25,ii,j,i)
     Gday_sun_high(ii,j,i) = Gday_sun_sort(ii_75,ii,j,i)
     Gday_diff_low(ii,j,i) = Gday_diff_sort(ii_25,ii,j,i)
     Gday_diff_high(ii,j,i) = Gday_diff_sort(ii_75,ii,j,i)
     Gday_zero_low(ii,j,i) = Gday_zero_sort(ii_25,ii,j,i)
     Gday_zero_high(ii,j,i) = Gday_zero_sort(ii_75,ii,j,i)
     TWday_low(ii,j,i) = TWday_sort(ii_25,ii,j,i)
     TWday_high(ii,j,i) = TWday_sort(ii_75,ii,j,i)
     Gtwday_low(ii,j,i) = Gtwday_sort(ii_25,ii,j,i)
     Gtwday_high(ii,j,i) = Gtwday_sort(ii_75,ii,j,i)

     Gday0d_sun_low(ii,j,i) = Gday0d_sun_sort(ii_25,ii,j,i)
     Gday0d_sun_high(ii,j,i) = Gday0d_sun_sort(ii_75,ii,j,i)
     Gday0d_diff_low(ii,j,i) = Gday0d_diff_sort(ii_25,ii,j,i)
     Gday0d_diff_high(ii,j,i) = Gday0d_diff_sort(ii_75,ii,j,i)
     Gday0d_zero_low(ii,j,i) = Gday0d_zero_sort(ii_25,ii,j,i)
     Gday0d_zero_high(ii,j,i) = Gday0d_zero_sort(ii_75,ii,j,i)
     TWday35d_low(ii,j,i) = TWday35d_sort(ii_25,ii,j,i)
     TWday35d_high(ii,j,i) = TWday35d_sort(ii_75,ii,j,i)

     end do
    end do
   end do


  ;==global average Gday and TWday at the hottest grid cells ranked by ensemble median
   wgt1d := ndtooned(wgt)
   opt = True
   opt@PrintStat = True
   ;print("summary of wgt1d"); 62% ocean grid cells, 38% land grids (14894 Mha)
   ;cstat := stat_dispersion(wgt1d, opt) ;detailed statistics

   ngrid = nlat*nlon ; global total grids 64800 of 180*360 (including ocean and _FillValue)
   lgrid = num(wgt1d.gt.0); 23968 land grids
   print("number of land grids: "+lgrid)
   print("number of global grids: "+ngrid)

   ;=if dim_pqsort_n sort x in decreasing order (top values on the left)
   gg_99 = tointeger(0.01*lgrid); hottest 1% land cells
   gg_999 = tointeger(0.001*lgrid)  ;grid index of 99.9th percentile (top 0.1%) land points
   gg_95 = tointeger(0.05*lgrid)  ;grid index of 95th percentile (top 5% land)
  
   ;--ensemble median at each grid cell 
   TWday2d := onedtond(ndtooned(TWday_med),(/nlev,ngrid/))
   Gtw2d := onedtond(ndtooned(Gtwday_med),(/nlev,ngrid/))
   hc2d := onedtond(ndtooned(hc_med),(/nlev,ngrid/))
   p2d := onedtond(ndtooned(P_med),(/nlev,ngrid/))

   TWday2d@_FillValue = -9.96921e+36 ; set a negative missing value so that they rank at the lower end
   Gtw2d@_FillValue = -9.96921e+36
   TWday2d_sort := TWday2d ; keep a copy of original daytime TW
   ip_TWday := dim_pqsort_n(TWday2d_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!
   TWday99g := TWday2d_sort(:,0:gg_99) ;
   TWday999g := TWday2d_sort(:,0:gg_999) 
   TWday95g := TWday2d_sort(:,0:gg_95)

   ;use permutation index of TWday to sort Gtw
   Gtw2d_sort := new((/nlev,ngrid/),"float")
   hc2d_sort := new((/nlev,ngrid/),"float")
   p2d_sort := new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      Gtw2d_sort(ii,:) = Gtw2d(ii,ip_TWday(ii,:))
      hc2d_sort(ii,:) = hc2d(ii,ip_TWday(ii,:))
      p2d_sort(ii,:) = p2d(ii,ip_TWday(ii,:))
   end do
   ;do not sort Gtw by itself, only need to find corresponding Gtw for each TWday
   Gtw99g := Gtw2d_sort(:,0:gg_99) ;global percentile of Gtw corresponding to the percentile of TWday
   Gtw999g := Gtw2d_sort(:,0:gg_999)
   Gtw95g := Gtw2d_sort(:,0:gg_95)
   hc99g := hc2d_sort(:,0:gg_99) ;global percentile of hc corresponding to the percentile of TWday
   hc999g := hc2d_sort(:,0:gg_999)
   hc95g := hc2d_sort(:,0:gg_95)
   P99g := p2d_sort(:,0:gg_99) ;global percentile of pa corresponding to the percentile of TWday
   P999g := p2d_sort(:,0:gg_999)
   P95g := p2d_sort(:,0:gg_95)
   
   Gday2d := onedtond(ndtooned(Gday_sun_med),(/nlev,ngrid/))
   Gday2d@_FillValue = -9.96921e+36 
   Gday2d_sort := Gday2d; missig values are sorted at the lower end
   ip_Gday_sun := dim_pqsort_n(Gday2d_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!
   Gday99g := Gday2d_sort(:,0:gg_99) ;99% percentile of Gday across the globe (land)
   Gday999g := Gday2d_sort(:,0:gg_999)
   Gday95g := Gday2d_sort(:,0:gg_95)

   Gday99g_sun_med = dim_avg_n(Gday99g,1) ; average Gday of hottest 1% grid cells
   Gday95g_sun_med = dim_avg_n(Gday95g,1)
   Gday999g_sun_med = dim_avg_n(Gday999g,1)
   TWday99g_med = dim_avg_n(TWday99g,1)
   TWday95g_med = dim_avg_n(TWday95g,1)
   TWday999g_med = dim_avg_n(TWday999g,1)
   Gtw99g_med = dim_avg_n(Gtw99g,1)
   Gtw95g_med = dim_avg_n(Gtw95g,1)
   Gtw999g_med = dim_avg_n(Gtw999g,1)
   hc99g_med = dim_avg_n(hc99g,1)
   hc95g_med = dim_avg_n(hc95g,1)
   hc999g_med = dim_avg_n(hc999g,1)
   P99g_med = dim_avg_n(P99g,1)
   P95g_med = dim_avg_n(P95g,1)
   P999g_med = dim_avg_n(P999g,1)

   ;--summary of low CI/percentile at each grid cell
   TWday2d := onedtond(ndtooned(TWday_low),(/nlev,ngrid/))
   Gtw2d := onedtond(ndtooned(Gtwday_low),(/nlev,ngrid/))

   ;use permutation index of TWday to sort Gtw and their low and high percentiles
   TWday2d_sort := new((/nlev,ngrid/),"float")
   Gtw2d_sort := new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      TWday2d_sort(ii,:) = TWday2d(ii,ip_TWday(ii,:))
      Gtw2d_sort(ii,:) = Gtw2d(ii,ip_TWday(ii,:))
   end do
   ;do not sort Gtw by itself, only need to find corresponding Gtw for each TWday
   TWday99g := TWday2d_sort(:,0:gg_99) ;
   TWday999g := TWday2d_sort(:,0:gg_999)
   TWday95g := TWday2d_sort(:,0:gg_95)
   Gtw99g := Gtw2d_sort(:,0:gg_99) ;global percentile of Gtw corresponding to the percentile of TWday
   Gtw999g := Gtw2d_sort(:,0:gg_999)
   Gtw95g := Gtw2d_sort(:,0:gg_95)

   Gday2d := onedtond(ndtooned(Gday_sun_low),(/nlev,ngrid/))
   ;Gday2d@_FillValue = -9.96921e+36
   ;Gday2d_sort := Gday2d; missig values are sorted at the lower end
   ;ip_Gday := dim_pqsort_n(Gday2d_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!
   Gday2d_sort := new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      Gday2d_sort(ii,:) = Gday2d(ii,ip_Gday_sun(ii,:))
   end do
   Gday99g := Gday2d_sort(:,0:gg_99) ;99% percentile of Gday across the globe (land)
   Gday999g := Gday2d_sort(:,0:gg_999)
   Gday95g := Gday2d_sort(:,0:gg_95)

   Gday99g_sun_low = dim_avg_n(Gday99g,1) ; average Gday of hottest 1% grid cells
   Gday95g_sun_low = dim_avg_n(Gday95g,1)
   Gday999g_sun_low = dim_avg_n(Gday999g,1)
   TWday99g_low = dim_avg_n(TWday99g,1)
   TWday95g_low = dim_avg_n(TWday95g,1)
   TWday999g_low = dim_avg_n(TWday999g,1)
   Gtw99g_low = dim_avg_n(Gtw99g,1)
   Gtw95g_low = dim_avg_n(Gtw95g,1)
   Gtw999g_low = dim_avg_n(Gtw999g,1)

   ;--summary of high CI/percentile at each grid cell
   TWday2d := onedtond(ndtooned(TWday_high),(/nlev,ngrid/))
   Gtw2d := onedtond(ndtooned(Gtwday_high),(/nlev,ngrid/))

   TWday2d_sort := new((/nlev,ngrid/),"float")
   Gtw2d_sort := new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      TWday2d_sort(ii,:) = TWday2d(ii,ip_TWday(ii,:))
      Gtw2d_sort(ii,:) = Gtw2d(ii,ip_TWday(ii,:))
   end do
   ;do not sort Gtw by itself, only need to find corresponding Gtw for each TWday
   TWday99g := TWday2d_sort(:,0:gg_99) ;
   TWday999g := TWday2d_sort(:,0:gg_999)
   TWday95g := TWday2d_sort(:,0:gg_95)
   Gtw99g := Gtw2d_sort(:,0:gg_99) ;global percentile of Gtw corresponding to the percentile of TWday
   Gtw999g := Gtw2d_sort(:,0:gg_999)
   Gtw95g := Gtw2d_sort(:,0:gg_95)

   Gday2d := onedtond(ndtooned(Gday_sun_high),(/nlev,ngrid/))
   ;Gday2d@_FillValue = -9.96921e+36
   ;Gday2d_sort := Gday2d; missig values are sorted at the higher end
   ;ip_Gday := dim_pqsort_n(Gday2d_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!
   Gday2d_sort := new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      Gday2d_sort(ii,:) = Gday2d(ii,ip_Gday_sun(ii,:))
   end do
   Gday99g := Gday2d_sort(:,0:gg_99) ;99% percentile of Gday across the globe (land)
   Gday999g := Gday2d_sort(:,0:gg_999)
   Gday95g := Gday2d_sort(:,0:gg_95)

   Gday99g_sun_high = dim_avg_n(Gday99g,1) ; average Gday of hottest 1% grid cells
   Gday95g_sun_high = dim_avg_n(Gday95g,1)
   Gday999g_sun_high = dim_avg_n(Gday999g,1)
   TWday99g_high = dim_avg_n(TWday99g,1)
   TWday95g_high = dim_avg_n(TWday95g,1)
   TWday999g_high = dim_avg_n(TWday999g,1)
   Gtw99g_high = dim_avg_n(Gtw99g,1)
   Gtw95g_high = dim_avg_n(Gtw95g,1)
   Gtw999g_high = dim_avg_n(Gtw999g,1)

   ;==diffuse scenario
   Gday2d := onedtond(ndtooned(Gday_diff_med),(/nlev,ngrid/))
   Gday2d@_FillValue = -9.96921e+36
   Gday2d_sort := Gday2d; missig values are sorted at the lower end
   ip_Gday_diff := dim_pqsort_n(Gday2d_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!
   Gday99g := Gday2d_sort(:,0:gg_99) ;grid cells >99th percentile (hottest 1% land grid cells)
   Gday95g := Gday2d_sort(:,0:gg_95)
   Gday999g := Gday2d_sort(:,0:gg_999)
   Gday99g_diff_med = dim_avg_n(Gday99g,1)
   Gday95g_diff_med = dim_avg_n(Gday95g,1)
   Gday999g_diff_med = dim_avg_n(Gday999g,1)

   Gday2d := onedtond(ndtooned(Gday_diff_low),(/nlev,ngrid/))
   ;Gday2d@_FillValue = -9.96921e+36
   ;Gday2d_sort := Gday2d; missig values are sorted at the lower end
   ;ip_Gday := dim_pqsort_n(Gday2d_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!
   Gday2d_sort := new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      Gday2d_sort(ii,:) = Gday2d(ii,ip_Gday_diff(ii,:))
   end do
   Gday99g := Gday2d_sort(:,0:gg_99) ;grid cells >99th percentile (hottest 1% land grid cells)
   Gday95g := Gday2d_sort(:,0:gg_95)
   Gday999g := Gday2d_sort(:,0:gg_999)
   Gday99g_diff_low = dim_avg_n(Gday99g,1)
   Gday95g_diff_low = dim_avg_n(Gday95g,1)
   Gday999g_diff_low = dim_avg_n(Gday999g,1)

   Gday2d := onedtond(ndtooned(Gday_diff_high),(/nlev,ngrid/))
   ;Gday2d@_FillValue = -9.96921e+36
   ;Gday2d_sort := Gday2d; missig values are sorted at the lower end
   ;ip_Gday := dim_pqsort_n(Gday2d_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!
   Gday2d_sort := new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      Gday2d_sort(ii,:) = Gday2d(ii,ip_Gday_diff(ii,:))
   end do
   Gday99g := Gday2d_sort(:,0:gg_99) ;grid cells >99th percentile (hottest 1% land grid cells)
   Gday95g := Gday2d_sort(:,0:gg_95)
   Gday999g := Gday2d_sort(:,0:gg_999)
   Gday99g_diff_high = dim_avg_n(Gday99g,1)
   Gday95g_diff_high = dim_avg_n(Gday95g,1)
   Gday999g_diff_high = dim_avg_n(Gday999g,1)


   ;==zero solar 
   Gday2d := onedtond(ndtooned(Gday_zero_med),(/nlev,ngrid/))
   Gday2d@_FillValue = -9.96921e+36
   Gday2d_sort := Gday2d; missig values are sorted at the lower end
   ip_Gday_zero := dim_pqsort_n(Gday2d_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!
   Gday99g := Gday2d_sort(:,0:gg_99) ;grid cells >99th percentile (hottest 1% land grid cells)
   Gday95g := Gday2d_sort(:,0:gg_95)
   Gday999g := Gday2d_sort(:,0:gg_999)
   Gday99g_zero_med = dim_avg_n(Gday99g,1)
   Gday95g_zero_med = dim_avg_n(Gday95g,1)
   Gday999g_zero_med = dim_avg_n(Gday999g,1)

   Gday2d := onedtond(ndtooned(Gday_zero_low),(/nlev,ngrid/))
   ;Gday2d@_FillValue = -9.96921e+36
   ;Gday2d_sort := Gday2d; missig values are sorted at the lower end
   ;ip_Gday := dim_pqsort_n(Gday2d_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!
   Gday2d_sort := new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      Gday2d_sort(ii,:) = Gday2d(ii,ip_Gday_zero(ii,:))
   end do
   Gday99g := Gday2d_sort(:,0:gg_99) ;grid cells >99th percentile (hottest 1% land grid cells)
   Gday95g := Gday2d_sort(:,0:gg_95)
   Gday999g := Gday2d_sort(:,0:gg_999)
   Gday99g_zero_low = dim_avg_n(Gday99g,1)
   Gday95g_zero_low = dim_avg_n(Gday95g,1)
   Gday999g_zero_low = dim_avg_n(Gday999g,1)

   Gday2d := onedtond(ndtooned(Gday_zero_high),(/nlev,ngrid/))
   ;Gday2d@_FillValue = -9.96921e+36
   ;Gday2d_sort := Gday2d; missig values are sorted at the lower end
   ;ip_Gday := dim_pqsort_n(Gday2d_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!
   Gday2d_sort := new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      Gday2d_sort(ii,:) = Gday2d(ii,ip_Gday_zero(ii,:))
   end do
   Gday99g := Gday2d_sort(:,0:gg_99) ;grid cells >99th percentile (hottest 1% land grid cells)
   Gday95g := Gday2d_sort(:,0:gg_95)
   Gday999g := Gday2d_sort(:,0:gg_999)
   Gday99g_zero_high = dim_avg_n(Gday99g,1)
   Gday95g_zero_high = dim_avg_n(Gday95g,1)
   Gday999g_zero_high = dim_avg_n(Gday999g,1)


   ;---average annual number of days with Gday>0 of above hottest 1% grid cells ranked by Gday
   ;ensemble median
   G0_avgdays_2d := onedtond(ndtooned(Gday0d_sun_med),(/nlev,ngrid/))
   G0_avgdays_2d@_FillValue = -9.96921e+36
   G0_avgdays_sort := G0_avgdays_2d; missig values are sorted at the lower end
   ip_G0day_sun := dim_pqsort_n(G0_avgdays_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!

   ;use permutation index of Gday_sun to sort
   ;G0_avgdays_sort = new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      ;G0_avgdays_sort(ii,:) = G0_avgdays_2d(ii,ip_Gday_sun(ii,:))
      wgt_sort = wgt1d(ip_G0day_sun(ii,:))
   end do
   G0_avgdays_99g := G0_avgdays_sort(:,0:gg_99)
   wgt_99g := wgt_sort(0:gg_99)
   wgt_99g := conform_dims(dimsizes(G0_avgdays_99g),wgt_99g,(/1/)) ;  [nlev,n99]
   G0_avgdays_99g_sun_med = tofloat(dim_sum_n(G0_avgdays_99g*wgt_99g, 1)/dim_sum_n(wgt_99g, 1)); area-weighted median

   G0_avgdays_2d := onedtond(ndtooned(Gday0d_diff_med),(/nlev,ngrid/))
   G0_avgdays_2d@_FillValue = -9.96921e+36
   G0_avgdays_sort := G0_avgdays_2d; missig values are sorted at the lower end
   ip_G0day_diff := dim_pqsort_n(G0_avgdays_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!
   ;use permutation index of Gday_diff to sort
   do ii=0,nlev-1
      ;G0_avgdays_sort(ii,:) = G0_avgdays_2d(ii,ip_Gday_diff(ii,:))
      wgt_sort = wgt1d(ip_G0day_diff(ii,:))
   end do
   G0_avgdays_99g := G0_avgdays_sort(:,0:gg_99)
   wgt_99g := wgt_sort(0:gg_99)
   wgt_99g := conform_dims(dimsizes(G0_avgdays_99g),wgt_99g,(/1/)) ;  [nlev,n99]
   G0_avgdays_99g_diff_med = tofloat(dim_sum_n(G0_avgdays_99g*wgt_99g, 1)/dim_sum_n(wgt_99g, 1))

   G0_avgdays_2d := onedtond(ndtooned(Gday0d_zero_med),(/nlev,ngrid/))
   G0_avgdays_2d@_FillValue = -9.96921e+36
   G0_avgdays_sort := G0_avgdays_2d; missig values are sorted at the lower end
   ip_G0day_zero := dim_pqsort_n(G0_avgdays_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!
   ;use permutation index of Gday_zero to sort
   do ii=0,nlev-1
      ;G0_avgdays_sort(ii,:) = G0_avgdays_2d(ii,ip_Gday_zero(ii,:))
      wgt_sort = wgt1d(ip_G0day_zero(ii,:))
   end do
   G0_avgdays_99g := G0_avgdays_sort(:,0:gg_99)
   wgt_99g := wgt_sort(0:gg_99)
   wgt_99g := conform_dims(dimsizes(G0_avgdays_99g),wgt_99g,(/1/)) ;  [nlev,n99]
   G0_avgdays_99g_zero_med = tofloat(dim_sum_n(G0_avgdays_99g*wgt_99g, 1)/dim_sum_n(wgt_99g, 1))

   TW35_avgdays_2d := onedtond(ndtooned(TWday35d_med),(/nlev,ngrid/))
   TW35_avgdays_2d@_FillValue = -9.96921e+36
   TW35_avgdays_sort := TW35_avgdays_2d; missig values are sorted at the lower end
   ip_TW35day := dim_pqsort_n(TW35_avgdays_sort, -2, 1) ;sort on dim=1 "grid" in decreasing order (-2)!
   ;use permutation index of TWday to sort
   ;TW35_avgdays_sort := new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      ;TW35_avgdays_sort(ii,:) = TW35_avgdays_2d(ii,ip_TWday(ii,:))
      wgt_sort = wgt1d(ip_TW35day(ii,:))
   end do
   TW35_avgdays_99g := TW35_avgdays_sort(:,0:gg_99)
   wgt_99g := wgt_sort(0:gg_99)
   wgt_99g := conform_dims(dimsizes(G0_avgdays_99g),wgt_99g,(/1/)) ;  [nlev,n99]
   TW35_avgdays_99g_med = tofloat(dim_sum_n(TW35_avgdays_99g*wgt_99g, 1)/dim_sum_n(wgt_99g, 1))

   ;=low 
   G0_avgdays_2d := onedtond(ndtooned(Gday0d_sun_low),(/nlev,ngrid/))
   ;use permutation index of Gday_sun to sort
   G0_avgdays_sort = new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      G0_avgdays_sort(ii,:) = G0_avgdays_2d(ii,ip_G0day_sun(ii,:))
      wgt_sort = wgt1d(ip_G0day_sun(ii,:))
   end do
   G0_avgdays_99g := G0_avgdays_sort(:,0:gg_99)
   wgt_99g := wgt_sort(0:gg_99)
   wgt_99g := conform_dims(dimsizes(G0_avgdays_99g),wgt_99g,(/1/)) ;  [nlev,n99]
   G0_avgdays_99g_sun_low = tofloat(dim_sum_n(G0_avgdays_99g*wgt_99g, 1)/dim_sum_n(wgt_99g, 1)); area-weighted lowian

   G0_avgdays_2d := onedtond(ndtooned(Gday0d_diff_low),(/nlev,ngrid/))
   ;use permutation index of Gday_diff to sort
   do ii=0,nlev-1
      G0_avgdays_sort(ii,:) = G0_avgdays_2d(ii,ip_G0day_diff(ii,:))
      wgt_sort = wgt1d(ip_G0day_diff(ii,:))
   end do
   G0_avgdays_99g := G0_avgdays_sort(:,0:gg_99)
   wgt_99g := wgt_sort(0:gg_99)
   wgt_99g := conform_dims(dimsizes(G0_avgdays_99g),wgt_99g,(/1/)) ;  [nlev,n99]
   G0_avgdays_99g_diff_low = tofloat(dim_sum_n(G0_avgdays_99g*wgt_99g, 1)/dim_sum_n(wgt_99g, 1))

   G0_avgdays_2d := onedtond(ndtooned(Gday0d_zero_low),(/nlev,ngrid/))
   ;use permutation index of Gday_zero to sort
   do ii=0,nlev-1
      G0_avgdays_sort(ii,:) = G0_avgdays_2d(ii,ip_G0day_zero(ii,:))
      wgt_sort = wgt1d(ip_G0day_zero(ii,:))
   end do
   G0_avgdays_99g := G0_avgdays_sort(:,0:gg_99)
   wgt_99g := wgt_sort(0:gg_99)
   wgt_99g := conform_dims(dimsizes(G0_avgdays_99g),wgt_99g,(/1/)) ;  [nlev,n99]
   G0_avgdays_99g_zero_low = tofloat(dim_sum_n(G0_avgdays_99g*wgt_99g, 1)/dim_sum_n(wgt_99g, 1))

   TW35_avgdays_2d := onedtond(ndtooned(TWday35d_low),(/nlev,ngrid/))
   ;use permutation index of TWday to sort
   TW35_avgdays_sort := new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      TW35_avgdays_sort(ii,:) = TW35_avgdays_2d(ii,ip_TW35day(ii,:))
      wgt_sort = wgt1d(ip_TW35day(ii,:))
   end do
   TW35_avgdays_99g := TW35_avgdays_sort(:,0:gg_99)
   wgt_99g := wgt_sort(0:gg_99)
   wgt_99g := conform_dims(dimsizes(G0_avgdays_99g),wgt_99g,(/1/)) ;  [nlev,n99]
   TW35_avgdays_99g_low = tofloat(dim_sum_n(TW35_avgdays_99g*wgt_99g, 1)/dim_sum_n(wgt_99g, 1))

   ;high
   G0_avgdays_2d := onedtond(ndtooned(Gday0d_sun_high),(/nlev,ngrid/))
   ;use permutation index of Gday_sun to sort
   G0_avgdays_sort = new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      G0_avgdays_sort(ii,:) = G0_avgdays_2d(ii,ip_G0day_sun(ii,:))
      wgt_sort = wgt1d(ip_G0day_sun(ii,:))
   end do
   G0_avgdays_99g := G0_avgdays_sort(:,0:gg_99)
   wgt_99g := wgt_sort(0:gg_99)
   wgt_99g := conform_dims(dimsizes(G0_avgdays_99g),wgt_99g,(/1/)) ;  [nlev,n99]
   G0_avgdays_99g_sun_high = tofloat(dim_sum_n(G0_avgdays_99g*wgt_99g, 1)/dim_sum_n(wgt_99g, 1)); area-weighted highian

   G0_avgdays_2d := onedtond(ndtooned(Gday0d_diff_high),(/nlev,ngrid/))
   ;use permutation index of Gday_diff to sort
   do ii=0,nlev-1
      G0_avgdays_sort(ii,:) = G0_avgdays_2d(ii,ip_G0day_diff(ii,:))
      wgt_sort = wgt1d(ip_G0day_diff(ii,:))
   end do
   G0_avgdays_99g := G0_avgdays_sort(:,0:gg_99)
   wgt_99g := wgt_sort(0:gg_99)
   wgt_99g := conform_dims(dimsizes(G0_avgdays_99g),wgt_99g,(/1/)) ;  [nlev,n99]
   G0_avgdays_99g_diff_high = tofloat(dim_sum_n(G0_avgdays_99g*wgt_99g, 1)/dim_sum_n(wgt_99g, 1))

   G0_avgdays_2d := onedtond(ndtooned(Gday0d_zero_high),(/nlev,ngrid/))
   ;use permutation index of Gday_zero to sort
   do ii=0,nlev-1
      G0_avgdays_sort(ii,:) = G0_avgdays_2d(ii,ip_G0day_zero(ii,:))
      wgt_sort = wgt1d(ip_G0day_zero(ii,:))
   end do
   G0_avgdays_99g := G0_avgdays_sort(:,0:gg_99)
   wgt_99g := wgt_sort(0:gg_99)
   wgt_99g := conform_dims(dimsizes(G0_avgdays_99g),wgt_99g,(/1/)) ;  [nlev,n99]
   G0_avgdays_99g_zero_high = tofloat(dim_sum_n(G0_avgdays_99g*wgt_99g, 1)/dim_sum_n(wgt_99g, 1))

   TW35_avgdays_2d := onedtond(ndtooned(TWday35d_high),(/nlev,ngrid/))
   ;use permutation index of TWday to sort
   TW35_avgdays_sort := new((/nlev,ngrid/),"float")
   do ii=0,nlev-1
      TW35_avgdays_sort(ii,:) = TW35_avgdays_2d(ii,ip_TW35day(ii,:))
      wgt_sort = wgt1d(ip_TW35day(ii,:))
   end do
   TW35_avgdays_99g := TW35_avgdays_sort(:,0:gg_99)
   wgt_99g := wgt_sort(0:gg_99)
   wgt_99g := conform_dims(dimsizes(G0_avgdays_99g),wgt_99g,(/1/)) ;  [nlev,n99]
   TW35_avgdays_99g_high = tofloat(dim_sum_n(TW35_avgdays_99g*wgt_99g, 1)/dim_sum_n(wgt_99g, 1))




  ;--------save ouput to netcdf and csv files


   lev1     =  levels
   lev1!0    = "level" ;named variables
   lev1@long_name = "global warming levels 1-4 degree Celsius"
  
   Gday_sun_med!0="lev"
   Gday_sun_med&lev = lev1
   Gday_sun_med!1="lat"
   Gday_sun_med&lat = lat
   Gday_sun_med!2="lon"
   Gday_sun_med&lon = lon
   copy_VarCoords(Gday_sun_med, Gday_diff_med)
   copy_VarCoords(Gday_sun_med, Gday_zero_med)
   copy_VarCoords(Gday_sun_med, Gday_sun_low)
   copy_VarCoords(Gday_sun_med, Gday_diff_low)
   copy_VarCoords(Gday_sun_med, Gday_zero_low)
   copy_VarCoords(Gday_sun_med, Gday_sun_high)
   copy_VarCoords(Gday_sun_med, Gday_diff_high)
   copy_VarCoords(Gday_sun_med, Gday_zero_high)
   copy_VarCoords(Gday_sun_med, TWday_med)
   copy_VarCoords(Gday_sun_med, TWday_low)
   copy_VarCoords(Gday_sun_med, TWday_high)
   copy_VarCoords(Gday_sun_med, Gtwday_med)
   copy_VarCoords(Gday_sun_med, Gtwday_low)
   copy_VarCoords(Gday_sun_med, Gtwday_high)
   copy_VarCoords(Gday_sun_med, hc_med)
   copy_VarCoords(Gday_sun_med, P_med)
   copy_VarCoords(Gday_sun_med, Taday_med)

   copy_VarCoords(Gday_sun_med, Gday0d_sun_med)
   copy_VarCoords(Gday_sun_med, Gday0d_diff_med)
   copy_VarCoords(Gday_sun_med, Gday0d_zero_med)
   copy_VarCoords(Gday_sun_med, Gday0d_sun_low)
   copy_VarCoords(Gday_sun_med, Gday0d_diff_low)
   copy_VarCoords(Gday_sun_med, Gday0d_zero_low)
   copy_VarCoords(Gday_sun_med, Gday0d_sun_high)
   copy_VarCoords(Gday_sun_med, Gday0d_diff_high)
   copy_VarCoords(Gday_sun_med, Gday0d_zero_high)
   copy_VarCoords(Gday_sun_med, TWday35d_med)
   copy_VarCoords(Gday_sun_med, TWday35d_low)
   copy_VarCoords(Gday_sun_med, TWday35d_high)


  ;===============save 3d data to netcdf
  dirNc  = dirPath2       ; where nc file will be saved
  filNc  = "Gday-gridcell-warminglevel-10yravg-ens-median-vcbc-outdoor.nc"

  pthNc  = dirNc+filNc
  system("/bin/rm -f "+pthNc)       ; remove any pre-existing file
  ncdf = addfile(pthNc ,"c")        ; open output netCDF file

  ncdf->Gday_sun_med  = Gday_sun_med
  ncdf->Gday_diff_med  = Gday_diff_med
  ncdf->Gday_zero_med  = Gday_zero_med
  ncdf->Gday_sun_low  = Gday_sun_low
  ncdf->Gday_diff_low  = Gday_diff_low
  ncdf->Gday_zero_low  = Gday_zero_low
  ncdf->Gday_sun_high  = Gday_sun_high
  ncdf->Gday_diff_high  = Gday_diff_high
  ncdf->Gday_zero_high  = Gday_zero_high
  ncdf->TWday_med  = TWday_med
  ncdf->TWday_low  = TWday_low
  ncdf->TWday_high  = TWday_high
  ncdf->Gtwday_med  = Gtwday_med
  ncdf->Gtwday_low  = Gtwday_low
  ncdf->Gtwday_high  = Gtwday_high
  ncdf->hc_med  = hc_med
  ncdf->P_med  = P_med
  ncdf->Taday_med  = Taday_med

  ncdf->Gday0d_sun_med  = Gday0d_sun_med
  ncdf->Gday0d_diff_med  = Gday0d_diff_med
  ncdf->Gday0d_zero_med  = Gday0d_zero_med
  ncdf->Gday0d_sun_low  = Gday0d_sun_low
  ncdf->Gday0d_diff_low  = Gday0d_diff_low
  ncdf->Gday0d_zero_low  = Gday0d_zero_low
  ncdf->Gday0d_sun_high  = Gday0d_sun_high
  ncdf->Gday0d_diff_high  = Gday0d_diff_high
  ncdf->Gday0d_zero_high  = Gday0d_zero_high
  ncdf->TWday35d_med  = TWday35d_med
  ncdf->TWday35d_low  = TWday35d_low
  ncdf->TWday35d_high  = TWday35d_high


  ;===================Save 1d data to csv

  ;spatial statistics and extent of impact
  csv_file1 = dirPath2+"Gday_hot1pct_sun_warminglevels_10yravg_ens_median-vcbc-outdoor.csv"
  system("rm -rf " + csv_file1)
  csv_file2 = dirPath2+"Gday_hot1pct_diff_warminglevels_10yravg_ens_median-vcbc-outdoor.csv"
  system("rm -rf " + csv_file2)
  csv_file3 = dirPath2+"Gday_hot1pct_zero_warminglevels_10yravg_ens_median-vcbc-outdoor.csv"
  system("rm -rf " + csv_file3)

  csv_file10 = dirPath2+"TWday_hot1pct_warminglevels_10yravg_ens_median-vcbc-outdoor.csv"
  system("rm -rf " + csv_file10)

  csv_file13 = dirPath2+"Gday0d_hot1pct_warminglevels_10yravg_ens_median-vcbc-outdoor.csv"
  system("rm -rf " + csv_file13)

;---Convert data array to 1D 
  ;1) spatial statistics of hottest 1% (Gday99g) or 5%(Gday95g) grid cells
  field_names1 := (/"Level", "Gday999g_sun_med", "Gday99g_sun_med", "Gday95g_sun_med", "Gday999g_sun_high", "Gday99g_sun_high", "Gday95g_sun_high", "Gday999g_sun_low", "Gday99g_sun_low", "Gday95g_sun_low"/)
  header1 := [/str_join(field_names1,",")/]
  write_table(csv_file1, "w", header1, "%s")
  alist1  := [/levels, Gday999g_sun_med, Gday99g_sun_med, Gday95g_sun_med, Gday999g_sun_high, Gday99g_sun_high, Gday95g_sun_high, Gday999g_sun_low, Gday99g_sun_low, Gday95g_sun_low /]
  format1 = "%g,%g,%g,%g,%g,%g,%g,%g,%g,%g"
  write_table(csv_file1, "a", alist1, format1)

  field_names1 := (/"Level", "TWday999g_med", "TWday99g_med", "TWday95g_med","Gtw999g_med","Gtw99g_med","Gtw95g_med", \ 
       "TWday999g_high", "TWday99g_high", "TWday95g_high","Gtw999g_high","Gtw99g_high","Gtw95g_high", \ 
       "TWday999g_low", "TWday99g_low", "TWday95g_low","Gtw999g_low","Gtw99g_low","Gtw95g_low","hc999g_med","hc99g_med","hc95g_med","P999g_med","P99g_med","P95g_med"/)
  header1 := [/str_join(field_names1,",")/]
  write_table(csv_file10, "w", header1, "%s")
  alist1  := [/levels, TWday999g_med, TWday99g_med, TWday95g_med, Gtw999g_med,Gtw99g_med, Gtw95g_med, \
       TWday999g_high, TWday99g_high, TWday95g_high, Gtw999g_high,Gtw99g_high, Gtw95g_high, \
       TWday999g_low, TWday99g_low, TWday95g_low, Gtw999g_low,Gtw99g_low, Gtw95g_low, hc999g_med,hc99g_med,hc95g_med,P999g_med,P99g_med,P95g_med /]
  format1 = "%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g"
  write_table(csv_file10, "a", alist1, format1)


  ;======diffuse scenario

  ;1) spatial statistics of hottest 1% (Gday99g) or 5%(Gday95g) grid cells
  field_names1 := (/"Level", "Gday999g_diff_med", "Gday99g_diff_med", "Gday95g_diff_med", "Gday999g_diff_high", "Gday99g_diff_high", "Gday95g_diff_high", "Gday999g_diff_low", "Gday99g_diff_low", "Gday95g_diff_low"/)

  header1 := [/str_join(field_names1,",")/]
  write_table(csv_file2, "w", header1, "%s")
  alist1  := [/levels, Gday999g_diff_med, Gday99g_diff_med, Gday95g_diff_med, Gday999g_diff_high, Gday99g_diff_high, Gday95g_diff_high, Gday999g_diff_low, Gday99g_diff_low, Gday95g_diff_low /]
  format1 = "%g,%g,%g,%g,%g,%g,%g,%g,%g,%g"
  write_table(csv_file2, "a", alist1, format1)


  ;======zero radiation scenario

  ;1) spatial statistics of hottest 1% (Gday99g) or 5%(Gday95g) grid cells
  field_names1 := (/"Level", "Gday999g_zero_med", "Gday99g_zero_med", "Gday95g_zero_med", "Gday999g_zero_high", "Gday99g_zero_high", "Gday95g_zero_high", "Gday999g_zero_low", "Gday99g_zero_low", "Gday95g_zero_low"/)

  header1 := [/str_join(field_names1,",")/]
  write_table(csv_file3, "w", header1, "%s")
  alist1  := [/levels, Gday999g_zero_med, Gday99g_zero_med, Gday95g_zero_med, Gday999g_zero_high, Gday99g_zero_high, Gday95g_zero_high, Gday999g_zero_low, Gday99g_zero_low, Gday95g_zero_low /]
  format1 = "%g,%g,%g,%g,%g,%g,%g,%g,%g,%g"
  write_table(csv_file3, "a", alist1, format1)



  ;=====annual number of days with Gday>0 at the hottest 1% grid cells
  field_names1 := (/"Level", "G0_avgdays_99g_sun_med", "G0_avgdays_99g_diff_med", "G0_avgdays_99g_zero_med", "TW35_avgdays_99g_med", "G0_avgdays_99g_sun_high", "G0_avgdays_99g_diff_high","G0_avgdays_99g_zero_high",  "TW35_avgdays_99g_high", "G0_avgdays_99g_sun_low", "G0_avgdays_99g_diff_low", "G0_avgdays_99g_zero_low", "TW35_avgdays_99g_low"/)
  header1 = [/str_join(field_names1,",")/]
  write_table(csv_file13, "w", header1, "%s")
  alist1  = [/levels, G0_avgdays_99g_sun_med, G0_avgdays_99g_diff_med, G0_avgdays_99g_zero_med, TW35_avgdays_99g_med, G0_avgdays_99g_sun_high, G0_avgdays_99g_diff_high, G0_avgdays_99g_zero_high, TW35_avgdays_99g_high, G0_avgdays_99g_sun_low, G0_avgdays_99g_diff_low, G0_avgdays_99g_zero_low,  TW35_avgdays_99g_low /]
  format1 = "%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g"
  write_table(csv_file13, "a", alist1, format1)


 end if; DATA






end 
